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Abstract. In the analysis of variable stars, the problem 
of sampling is central. This article focusses on the deter¬ 
mination of the Nyquist frequency. It is well defined in the 
case of regular sampling. However, the time series of vari¬ 
able stars observations are generally unevenly sampled. 
Fourier analysis using the spectral window furnishes some 
clues about the equivalent Nyquist frequency in the ir¬ 
regular case. Often it is pushed very high, and thus very 
short periods can be detected. A specific example is shown, 
drawn from MACHO databases. 
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1. Introduction 

When variable stars observations are analyzed, two ma¬ 
jor problems arise. First, there are constraints on the 
sampling times. They can’t always be chosen, because of 
weather conditions or other observing restrictions. Sec¬ 
ondly, the stars, even if they are periodic, can have an ex¬ 
tremely wide range of periods and behaviors. Therefore, 
astronomers would like to know up to which higher limit 
it makes sense to search for frequency. 

In the literature, there are contradictions regarding 
the Nyquist frequency for irregularly sampled data. Some 
authors (Press et al. 1992, Horne et al. 1986) identify it 
with l/2(5t where St is the “average” sampling rate. The 
Nyquist frequency is also often identified with v = l/2s, 
where s is the smallest time interval in the sample (Scar- 
gle 1982, Roberts et al. 1987). Nevertheless, Press et al. 
(1992) and Roberts et al. (1987) remarked that frequencies 
can be detected above their quoted values. 

For evenly sampled data, the Nyquist frequency is de¬ 
fined as l/26t, where St is the time sampling step. We will 
show how this quantity can be extended to the irregular 
situation. In most cases, the Nyquist frequency is much 
larger than usually thought. The present work is particu¬ 
larly relevant to large databases that are becoming avail¬ 
able (MACHO, OGLE, EROS, HIPPARCOS, etc...). 
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2. Nyquist frequency for irregular sampling 

Let /(t) be a signal function of time t. We suppose that 
the power spectrum of the function / is not variable in 
time (it should be noted that this hypothesis is stronger 
than stationarity). This is the case for periodic functions. 
Otherwise, the Nyquist frequency is only defined very lo¬ 
cally and loses its sense for a given irregular sampling. 
Moreover, a spectral analysis of an unstable signal which 
is undersampled may lead to wrong conclusions. 

If / is sampled at time ti we are left with {fi,ti) 
i = 1,..., A, where N is the total number of measure¬ 
ments. In the frequency domain. Deeming (1975, 1976) 
showed how the irregular discretization affects the power 
spectrum : the spectrum results from the con¬ 

volution between the real spectrum F{v) and the spec¬ 
tral window Gi<i{v)- is the Fourier transform of 

g{t) — '^S{t — ti). That is Fi\f{h') = F(:/) (g) Gn{v). 

We can state the result as follows. Let[]p be the largest 
value such that Vt^, ti = ti + riip , where G N. p is a 
kind of greatest common divisor (gcd) for all (ti—ti). The 
Nyquist Frequency then is: 


^Ny 


1 1 
2p ~ 2s 


( 1 ) 


In most practical cases p <C s, therefore, even with a 
strong but random under-sampling, very high frequencies 
can be detected. 


3. Proof 

The key function is the spectral window: 

= I r ( 2 , 

We have Gn{v) < 1 and Gn{v) = Ifov u = m/p, m G N. 
Replacing tk with n^p + ti, we notice that this function 
is periodic with period 1/p (each individual term in the 
sum has this same periodicity). Moreover as g{t) is real, 
the function Gfq{v) is symmetric around the origin. These 

^ Note that if the ti are irrational, p may not exist and should 
then be taken as zero. 
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properties of periodicity and symmetry imply that Gn{v) 
is also symmetric around {2k + l)/2p. For A; = 0, we have 
l/2p, which is just the Nyquist frequency. For the power 
spectrum, the proof is exactly the same. We can think of 
the set ti as a regular sampling with step p, where some 
(most) data have been dropped. The more random the 
data is, the smaller p will be and so the higher As 
in the regular case, the continuous signal f{t) should not 
contain any component above Otherwise, it will be 
mirrored into 2vt^y — v, and, without further assumptions, 
it will not be distinguishable from the mirror image. Of 
course, no frequency component filtered during the mea¬ 
surement process can be recovered by the random sam¬ 
pling. In particular, if Atexp is the exposure time of ev¬ 
ery measures, then in the frequency domain, the signal is 
filtered by sin^(7rAtexpJ^)/(7i'Atexp^^)^ and not much infor¬ 
mation will be recovered above t'max = l/(KAtexp), k = 
2 or 3. 
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4. Simulation 


The following numerical experiment illustrates the theory 
presented above. We chose the times as 


[1000-rndj , ^ 

Too 


,20 


( 3 ) 


where rnd are random numbers uniformly distributed 
on [0,1[. Then, we built up a signal as follows f{ti) = 
3sin(27rr'oii), with i^o = 43 (period = .0233...). 


Time 

Signal 

Time 

Signal 

0.00 

-2.69 

5.59 

0.88 

0.52 

0.70 

6.12 

-2.55 

1.07 

-2.77 

7.27 

2.91 

1.33 

-2.21 

7.84 

-2.86 

1.77 

-2.91 

8.13 

2.98 

2.00 

-2.69 

8.32 

1.14 

2.77 

-2.91 

8.49 

-2.99 

3.31 

0.14 

8.54 

-1.80 

4.30 

-1.40 

9.00 

-2.69 

4.49 

-2.99 

9.29 

2.39 


The smallest time interval is s = 0.05, so vq is well above 
l/2s, but p is 0.01. As expected, we observe that the 
Nyquist frequency is at 50, and that there is no strong 
aliasing peak before v = 100 (figure 1). 

If we had rounded the times tk to 0.001 in the above 
example and had p at 0.001, we would have found at 
500. 

5. Determination of p 

Theoretically, p could be found using an extension to more 
than 2 numbers of the Euclid gcd algorithm. But in prac¬ 
tice, the tk are not known to infinite precision, and the 
Euclid algorithm is numerically unstable for non integer 


Fig. 1. Spectral window Gn and power spectrum Fpf for 
the simulation described in the text. Both functions are 
clearly symmetrical around r'Ny = l/2p = 50. No peak 
appears for Gn{v) ai v = 1/s. 


numbers. A practical solution could then be to stop the al¬ 
gorithm when the smallest number is smaller than some e. 
Otherwise, the spectral window is probably the most ef¬ 
ficient way to find p in this fuzzy environment, p is the 
inverse of the smallest v for which Gn{i') is greater than 
some value just below one. 


6. Application 

Actually, most data are obtained at approximately equal 
time intervals. The fact that there is some randomness 
may raise substantially the Nyquist frequency. As a con¬ 
spicuous example, we can cite Minniti et al. (1998) who 
are able to find periods for 6 Scuti stars with MACHO ob¬ 
servations of bulge fields. Indeed, the data are separated 
by about one day and the periods found are of few hours. 
We give an example for the 5 Scuti star 162.25348.3066. 
The histogram of the time intervals is given in figure 2, 
its Fourier transform is shown in figure 3, and finally the 
folded curve is presented in figure 4, it leaves no doubts 
that the period is the right one. The detection limit is 
dominated by the exposure time (150 seconds), that is 
about JZmax = 192 days”^ if k = 3 (Pmin = 7-5 min). 

With the HIPPARCOS mission the smallest time in¬ 
terval is 20 minutes. A search for shorter periods was un¬ 
dertaken (Eyer, 1998). Unfortunately no results came out. 
We can invoke that often the spectral window is nearly 
symmetric around 1/40 min“^ (this is due to the fact that 
observations are often found in sequences, inducing strong 
aliasing peaks), or that short periods are rarely stable over 
three years (either period change or phase shift can occur). 
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dt [days] 

Fig. 2. Histogram of the time intervals. The smallest is 
around 0.8 days. In total, 122 measurements were avail¬ 
able. 



Fig. 3. Spectral window Gn and power spectrum Fjv of 
the S Scuti star 162.25348.3066 from MACHO. Although 
there are rather strong aliasing peaks, the frequency of 
1/0.0854259 days”^ is unambiguous. 


7. Reducing the amplitude of peaks in Gjv(^') 

The spectral window also helps to decide when new mea¬ 
surements should be made to reduce the annoying spikes 
in the spectral window. We make use of the simulated 
data in section 4 as an example. A practical procedure 
could be to determine the polar coordinates p, 9 (point 
(1) in figure 5, where we used as a concrete example 
the simulated data of section 4) of the complex number 
(1/A^) ^^^1 exp(z27r:/;,tfc) for the annoying frequency Vh- 
Point (1) is the center of gravity of the observing points. 
The value of is the distance of point (1) to the 

origin. Then, the new measurement should be taken at 


tN+l 



TT -I- 
27r ) 


1 


(4) 
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Fig. 4. The folded curve of the 6 Scuti star 162.25348.3066 
from MACHO (period = 0.0854259 days). 
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Fig. 5. The points exp{i2Tn/htk) are displayed in the 
complex plane (small open squares) for the frequency 
Vh = 3.99 corresponding to the highest peak of Gn{v) 
{1 < V < 50, in figure 1). The black square (1) represents 
the center of gravity of these points, while the open square 
(2) opposed to (1) on the circle, is the best position for 
tN+i (modulo 2Trvh) to reduce the amplitude of Gat at Vh- 


"a (2) 


( 1 ) 


-0.5 0 0,5 1 

cos(27ti/ t^) 


for some integer value k compatible with the other con¬ 
straints. This produces the largest reduction of the peak. 
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